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ABSTRACT 

In this paper we describe a Bayesian statistical method designed to infer the mag- 
netic properties of stars observed using high-resolution circular spectropolarimetry 
in the context of large surveys. This approach is well suited for analysing stars for 
which the stellar rotation period is not known, and therefore the rotational phases 
of the observations are ambiguous. The model assumes that the magnetic observa- 
tions correspond to a dipole oblique rotator, a situation commonly encountered in 
intermediate and high-mass stars. Using reasonable assumptions regarding the model 
parameter prior probability density distributions, the Bayesian algorithm determines 
the posterior probability densities corresponding to the surface magnetic field geom- 
etry and strength by performing a comparison between the observed and computed 
Stokes V profiles. 

Based on the results of numerical simulations, we conclude that this method yields 
a useful estimate of the surface dipole field strength based on a small number (i.e. 1 
or 2) of observations. On the other hand, the method provides only weak constraints 
on the dipole geometry. The odds ratio, a parameter computed by the algorithm that 
quantifies the relative appropriateness of the magnetic dipole model versus the non- 
magnetic model, provides a more sensitive diagnostic of the presence of weak magnetic 
signals embedded in noise than traditional techniques. 

To illustrate the application of the technique to real data, we analyse seven ES- 
PaDOnS and Narval observations of the early B-type magnetic star LP Ori. Insufficient 
information is available to determine the rotational period of the star and therefore 
the phase of the data; hence traditional modelling techniques fail to infer the dipole 
strength. In contrast, the Bayesian method allows a robust determination of the dipole 
polar strength, Ba = 911^111 G. 

Key words: stars: magnetic fields- stars: early- type - techniques: polarimetric - 
methods: statistical. 



1 INTRODUCTION 

The evolution of a massive star is strongly determined by its 
rotation, as well as the mass lost through its stellar wind, 
both of which can be strongly influenced by the presence of a 
magnetic field (e.g. Townsend et al.|2010 Wade et al.|2011| ). 
A magnetic field can furthermore couple different layers of 
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a star's interior, thereby modifying internal differential ro- 



tation and circulation currents (Maeder & Meynet 2005). 
If a field has a large-scale component that extends outside 
the stellar surface, it can also channel the outflowing stel- 
lar wind, creating a structured wind - a magnetosphere - 
which will modify the rate and geometry of mass loss along 
with its observational properties ( Landstreet Borra|1978 



Shore Bro^ fToOQ) [Babel Montmerle||1997a|b| ). Fur- 



thermore, if the held couples the rotating surface of the star 
with the outflowing stellar wind, both effects will result in 
angular momentum loss (via the outflowing stellar wind), 
which differs strongly from that of a non-magnetic star (ud- 
jDoula et aL|[2b09). As angular momentum and mass loss 
are cornerstone inputs in stellar evolution calculations, it is 
crucial that the effect of magnetic flelds in massive stars be 



2 V. Petit & G. A. Wade 



understood properly. For example, evolutionary tracks and 
isochrones can be used to interpret large datasets associ- 
ated with OB associations like the VLT-FLAMES surveys 
of massive stars ( Hunter et al.|2008t . 

Over the last decade, our knowledge of the properties 
of massive star magnetic fields has significantly improved, 
in large part due to a new generation of powerful high- 
resolution spectropolarimetric instrumentation. Tradition- 
ally, stellar magnetic fields were modelled using measure- 
ments of the longitudinal magnetic field of the star, yielding 
a single quantity: the strength of the disc-integrated line-of- 
sight component of the magnetic field (e.g. Bag nulo et aLj 
|2006|. In the absence of a field detection, the interpretation 



of such data was entirely dependent on the (unknown) stel- 
lar and magnetic geometry, and therefore highly ambiguous. 
However, modern high-resolution spectropolarimetry yields 
Doppler-broadened, velocity-resolved Stokes profiles mea- 
sured across spectral lines, which encode additional infor- 
mation about the field geometry. These data therefore allow 
a more robust inference of the field characteristics. 

Remarkably detailed information about the strength 
and local/global structure of stellar magnetic fields can be 
extracted from high-resolution, phase-resolved spectropo- 
larimetric datasets acquired in two or four Stokes parameters 
(i.e. Stokes IV or Stokes IVQU) using techniques such as 
Magnetic Doppler Imaging (e.g. [Kochukhov &: W ade 2010). 
However, detailed modelling of this sort relies on exten- 
sive high-quality datasets that demand significant telescope 
time, and which can only be obtained for a small number 
of stars. The lower-quality data that can be obtained for 
stars with less suitable observational properties can still be 
approximately investigated using parametric models; never- 
theless, even such a simple approach often requires approx- 
imately a dozen observations per star. 

Large observing programs, such as the Magnetism in 
Massive Stars (MiMeS) project, have dedicated significant 
resources to survey large samples of massive stars in search 
of magnetic fields ( Wade et al.|2011 ). Such surveys typically 
acquire a small number of Stokes V observations (typically 
1-3) per star, for a large number of stars. An outstanding 
problem concerns how to extract useful information about 
the magnetic field of a star from such a limited data set. 

In this paper we describe an approach to constrain the 
magnetic field strength of a star from small numbers of 
high-resolution Stokes V observations, assuming the dipole 
oblique rotator paradigm, and expressed using the formal- 
ism of Bayesian statistics. In Sect. [2] we describe briefly the 
model used to synthesise the emergent local circular po- 
larisation produced by the Zeeman effect under the weak- 
field approximation. We present the disk-integrated emer- 
gent Stokes V profiles obtained for a dipolar magnetic topol- 
ogy of a rotating star. In Sect. [3] we introduce the Bayesian 
algorithm used to compare a set of high-resolution Stokes 
V observations with a grid of synthetic line profiles. Sect. [4] 
presents the application of the Bayesian algorithm to sim- 
ulated data. We demonstrate that the Bayesian odds ratio 
can help detect a weak magnetic signal, by quantifying which 
target should be re-observed. It can also discriminate, with 
a few additional observations, whether an observation has a 
noise pattern that by chance looks like a magnetic signal ver- 
sus a real magnetic signal. We also show that the Bayesian 
algorithm provides a meaningful upper limit on the mag- 



netic strength in the case of a non-detection, and is able to 
reliably estimate the dipole field strength in the case of a 
detection. Some limited constraints can also be obtained for 
the geometry of the dipole. We compare the Bayesian di- 
agnostics with the traditional global longitudinal field and 
signal detection probability diagnostics. Finally, in Sect. [5] 
we present the Bayesian results for the magnetic B-type star 
LP Ori (P etit et al.^ 2008), obtained with observations of var- 
ious signal-to-noise ratios. 



2 LINE SYNTHESIS 
2.1 Local profiles 

Splitting of spectral lines due to the Zeeman effect corre- 
sponds to about lkms~^ per kG of surface field modulus. 
This implies that even relatively strong fields are difficult 
to detect reliably from Zeeman splitting in the spectrum 
of most stars. The situation is particularly challenging in 
the case of hot stars, where thermal broadening (of order a 
few kms~^), turbulent broadening (of order a few tens of 
kms~^) and rotational broadening (potentially a few hun- 
dreds of kms~^) combine to fully obscure any modification 
of the line profile due to the Zeeman effect. 

The Zeeman effect provides us with a second useful tool, 
in the form of the polarisation of the Zeeman components. 
Zeeman components exist in two different types: tt compo- 
nents, spread symmetrically about the zero- field wavelength 
of the spectral line and a components, with symmetric wave- 
length offsets to the red and blue of the zero-field wave- 
length. If the external magnetic field is oriented parallel to 
the observer's line of sight (a longitudinal field), the tt com- 
ponents vanish, while the two groups of a components have 
opposite circular polarisations. In a field aligned perpendicu- 
lar to the line of sight (a transverse field), the tt components 
are linearly polarised perpendicular to the field direction, 
while the a components are linearly polarised parallel to 
the field direction. Therefore, spectropolarimetric observa- 
tions are an extremely useful tool for detection and char- 
acterisation of both the strength and orientation of stellar 
magnetic fields. 

In order to correctly predict the polarised spectrum that 
will emerge from a stellar atmosphere, one must perform ra- 
diative transfer in an anisotropic medium, where the prop- 
agation properties depend on the propagation directions. In 
the weak- field regime, we can treat the anisotropic absorp- 
tion and dispersion profiles as perturbations of the isotropic 
case ( |Landi degl 'Innocent i Landolfi|2004| . In that case, at 
first order, only the longitudinal component of the magnetic 
field contributes to the polarisation, and therefore only cir- 
cular polarisation is predicted. The contribution of the trans- 
verse field, and the occurrence of linear polarisation, is a 
second-order effect. Therefore, circular polarisation (Stokes 
V spectra) signatures are generally stronger, and used for 
large surveys as such observations provide the best detec- 
tion threshold. 

Although there are some elaborate codes that are able 
to accurately synthesise the detailed circular polarisation 
profile of an arbitrary spectral transition (for a review see 
Wade et al.|2001 ), we will use the weak- field approximation, 
in order to draw a simple picture of the Stokes V spectrum 
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emerging from a star. The 1st order solution for the circular 
polarisation to the polarised radiative transfer equations is 
given by: 



V(v) 



^-C^effAo^ll 



dl{v) 
dv 



(1) 



The Stokes V profile V{v) emerging from a point at the sur- 
face of a star (referred to as a local profile) has the same 
shape as the derivative of the local intensity profile I{v), 
scaled by the longitudinal component of the magnetic field 
B\\ at that point, by the wavelength Ao of the transition, and 
by the effective Lande factor ^eff of the transition. The ef- 
fective Lande factor represents the separation of a triplet 
Zeeman pattern that approximates a more complex Zee- 
man pattern, and can be calculated under LS coupling, or 
taken from atomic experiments. Although we will be us- 
ing the weak field approximation throughout this paper, the 
Bayesian algorithm can use synthetic profiles computed with 
any spectrum synthesis code. 

A multi-line approach will usually be applied - in order 
to increase the signal to noise ratio (s/n) - to the data for 
which a Bayesian approach will be useful. Although these 
multi-line techniques have been developed and refined for 
more than two decades now (e.g. iBorr a et al1|1981| |Semel| 
& Li 1996| ), the most widely-used approach is the Least 
Squares Deconvolution (LSD) method introduced by Do- 
nati et aT] (fl997). Under the weak field approximation, the 
LSD method assumes that all the spectral lines have the 
same shape, and that this shape is scaled by the line depth 
for the intensity line profiles (Stokes /) and by the product 
of the line depth, the wavelength and the effective Lande 
factor for the circularly polarised line profile (Stokes V). 
The LSD method therefore already approximates the Zee- 
man pattern of a line by a triplet whose separation is given 
by an effective Lande factor. For a complete discussion on 
the circumstances under which a LSD profile can be treated 



as a single spectral line, see Kochukhov et al. (2010) 



2.2 Disk integration 

In order to predict the Stokes V spectrum of a given stellar 
spectral line, we must take into account the contribution 
from every point on the visible hemisphere. In this paper, 
we consider a single spectral line of arbitrary depth and a 
width corresponding to the sum of the spectral resolution 
and an ad hoc local broadening width. The contribution 
of each point on the visible hemisphere is dictated by the 
effective surface area and a wavelength-independent limb- 
darkening of the form: 



lc,0 



e + e cos ( 



(2) 



where Ic,o is the intensity at the stellar disk centre, 6 is 
the angle between the normal to the surface and the line of 
sight, and e is the limb-darkening factor (Gray 1992). We 



adopt e = 0.6 for the rest of this paper. The surface velocity 
field of the stellar disk is defined by the projected equatorial 
velocity vsini assuming rigid rotation. 

We consider a simple magnetic topology with a dipolar 
field of strength Bp, as the vast majority of intermediate 
mass and high mass stars' magnetic fields appear to be pre- 
dominantly dipolar (e.g. Borra Landstreet||1980 Auriere 
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Figure 1. Radial (top) and longitudinal (bottom) surface mag- 
netic field for a dipolar field seen pole-on (left) and equator-on 
(right). For the radial field, red, blue and green colours represent 
field oriented away from, into, and parallel to the star surface, re- 
spectively. For the longitudinal field, red, blue and green colours 
represent fields oriented toward, away from and perpendicular to 
the observer's line of sight, respectively. The black grid represent 
the rotation reference frame, with an inclination i = 90°. The 
obliquity of the field is ^ = 90° and the pole-on and equator-on 
configurations would correspond to rotational phases (/? = 0° and 
= 90°, respectively. 



et al.|200'7 ^. The orientation of the dipole with respect to the 
observer's reference frame is described by the inclination i of 
the rotational axis to the line of sight, the rotational phase 
(p and the obliquity /3 of the magnetic axis with respect to 
the rotational axis. 

Figure [l] shows the radial field (top) and the longitu- 
dinal field (bottom) for a dipole field seen positive pole-on 
(left) and magnetic equator-on (right). The radial field is 
defined by its orientation with respect to the stellar surface 
normal. For example, when we are looking at the positive 
pole (left), the radial field is at its maximum (red) in the 
middle of the stellar disk. On the edge of the disk, the field 
is parallel to the surface, and the radial field is null (green). 

The longitudinal field is the important quantity for the 
polarised radiative transfer. The longitudinal field represents 
the projection of the magnetic field vectors with respect 
to the observer's line of sight (perpendicular to the paper 
plane) . If we look at the magnetic field seen pole-on (bottom 
left), all of the magnetic vectors that are oriented toward the 
observer are colour-coded in red shades, green corresponds 
to a null longitudinal field component, and magnetic vectors 
oriented away from the observer are colour-coded in blue 
shades. As seen in Eq. ([T]), two longitudinal magnetic fields 
of same magnitude but opposite signs will produce inverse 
local Stokes V profiles. However, a good fraction of the field 
components oriented away from the observer are located on 
the hidden hemisphere of the star. There is therefore a net 
positive longitudinal component on the visible hemisphere 
(i.e. the visible hemisphere is more red than blue). In terms 
of the total emergent Stokes V, positive local profiles (lo- 
cal profiles for which the longitudinal field is positive) will 
dominate the total line profile. This effect is even more pro- 
nounced when limb darkening is taken into account, as the 
edges of the disk contribute less to the total brightness. 
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We now turn our attention to the case where we are 
looking at the magnetic equator (bottom right). In prin- 
ciple, there are as many magnetic vectors pointing toward 
the observer as away. The global longitudinal component is 
therefore null, as every magnetic vector on the left side of 
the stellar disk has its opposite on the right side of the disk. 
In terms of the emerging Stokes V profile, each positive lo- 
cal Stokes V profile will be cancelled out by the negative 
local Stokes V profile of same amplitude coming from the 
opposite side of the visible stellar disk, even when consid- 
ering limb-darkening. The resulting Stokes V profile should 
therefore be null when the global longitudinal component is 
null. This is effectively the case, if the star is not rotating. 

If the star is rotating around a rotation axis oriented to- 
ward the top of the paper (as represented by the black grid), 
the line of sight component of that rotation will introduce 
some Doppler shifts to different points on the stellar surface. 
For example, if the star rotates from left to right, the point 
situated to the extreme left of the stellar disk will be shifted 
by —vsini. The point situated at the extreme right of the 
stellar disk, whose local Stokes V profile would in principle 
cancel out the one produced by the leftmost point, will now 
be shifted by +vs\ni. Therefore, the rotation is able to sep- 
arate in the velocity space the local Stokes V line profiles 
that would otherwise cancel out, and a net circular polarisa- 
tion profile can be seen even if the global longitudinal field is 
null. It has been shown that a v sini as low as a few kms~^ 
is enough for instruments able to resolve this effect in the 
line profile CWad e et al.]|2000b| . 

Figure |2] shows an example of the shapes and relative 
amplitudes of Stokes V profiles emerging from a star rotat- 
ing at 50kms~^ (top). The radial and longitudinal fields 
are shown (middle), as well as the global longitudinal field 
curve (bottom). The phases indicated by red dots on the 
curve indicate the rotational phases corresponding to each 
shown profile. We can see that for this dipole configuration 
{i — 90°, (3 — 90°), the emerging Stokes V profile has an 
amplitude as strong when the global longitudinal field is 
null (2nd and 4th phases) than when it is at its maximum 
(1st, 3rd and last phases). We note the Stokes V profiles at 
= G are symmetric around the centre of the line. If we 
were to observe such profiles with an instrument that does 
not resolve the line profiles, we would indeed obtain a null 
global longitudinal field, even if the profiles are as clearly 
detectable with high-resolution observations than when the 
longitudinal field is at its maximum. 

Figure [3] shows a second dipole configuration (i — 45°, 
P = 45°), where the rotation is not able to separate com- 
ponents that will cancel out. When the longitudinal field is 
null (middle profile, we are looking at the magnetic equa- 
tor), we can see that as the star rotates from left to right, 
all the points situated on the left side of the stellar disk 
will be shifted to the blue. However, each half of the visible 
hemisphere contains as many field vectors oriented toward 
the observer as oriented away from the observer. Therefore, 
all the red-shaded points situated on the left side of the 
stellar disk will have a corresponding inverse profile that 
will be Doppler shifted by the same amount, and the can- 
cellation will occur. This particular symmetry occurs when 
P + i ^ 90°, at a phase where we are looking at the magnetic 
equator. 

Due to the intrinsic symmetry of a dipole field, as well 
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Figure 2. Top: Shape of the Stokes V line profiles for five rota- 
tional phases of a star with a dipolar field with i = (3 = 90°. Mid- 
dle: Corresponding radial field and longitudinal field components 
(colour schemes are as indicated in Figure The rotation axis 
inclination i = 90° is represented by the black grid and the mag- 
netic field pole is perpendicular to the rotation axis (/3 = 90°). 
Bottom: Global longitudinal field curve, normalised to the dipole 
field strength. The dots show the phase of the profiles. Note how 
the Stokes V profile do not disappear when the global longitudinal 
field goes to zero. 
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Figure 3. Same as Figure[2]for a rotation axis inclined hy i = 45° 
and a magnetic pole of obliquity (3 = 45°. Note how the Stokes V 
signal disappears when the global longitudinal field goes to zero 
because i + /3 ~ 90°. 



as the fact that circular polarisation is only sensitive to the 
longitudinal field components, some parameter degeneracy is 
encountered for the resulting Stokes V profiles. For example, 
if a dipolar field of parameter i, /3 and Lp produces a profile 
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/, the following symmetries (or anti-symmetries) occur: 
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(3) 



3 THE BAYESIAN ALGORITHM 

Bayesian statistics have proven to be an useful tool to find 



plausible models to explain astronomical data (e.g. Tuomi| 
et al.||2011||Arregui h Asensio RamospOllt . A probability 
density distribution p{Hi\I) provides a quantitative encod- 
ing of the plausibility of a certain hypothesis (Hi), given 
our current state of knowledge (/). The main pathway to 
Bayesian statistics |jaynes 2003) is the Bayes theorem: 



p{Hi\D,I) 



p{Hi\I)piD\Hi,I) 
V{D\I) 



(4) 



This theorem states that the posterior probability den- 
sity p{Hi\D,I) of an hypothesis, given a new dataset D 
and the current state of knowledge, is equal to the prod- 
uct of prior knowledge of the plausibility of the hypothe- 
sis p{Hi\I) and the likelihood p{D\Hi,I) of obtaining the 
new dataset if the hypothesis is true. The global likelihood, 
p{D\I) = '^^p(Hi\I)p{D\Hi, I), SiCts as a normalisation con- 
stant. The act of summing (or integrating in the case of a 
continuous hypothesis set) is called marginalisation. 

One difficulty of modelling the Stokes V spectra of stars 
observed as part of a survey like MiMeS is that most of 
the time the rotational phase of the star at the moment of 
the observation is not known, because the rotational pe- 
riod is generally unknown. Our approach is therefore to 
use Bayesian probability nomenclature to find which set of 
phase-independent configurations B = {Bp,i,l3) can repro- 
duce the observations. In other words, while the geometry 
B of the field must stay the same for each observation n of 
a given star, the phase cpn may in principle have any value. 

The hypothesis to be tested is the presence of an oblique 
dipolar magnetic field in a particular star. The predic- 
tions of this hypothesis are represented by the model Mi, 
parametrized with the parameters B {=[Bp, i, /S]) and $ 
{=[(Pi..(Pn]), the latter representing a set of phases associ- 
ated with a set of Stokes V observations V {=[Di..Dn])- 
The Bayes theorem tells us that the joint posterior proba- 
bility density for the parameters, assuming the veracity of 
the model Mi, is 



p(i3,$|P,Mi 



p{B,^\Mi)p{V\B,^,Mi) 
p{V\M^) 



(5) 



Our prior knowledge of the probability density for each 
model parameter, which can be as simple as its expected 
range, can be encoded in the prior term p{B^ $|Mi), whereas 
the information brought forward by the new observations is 
reflected in the likelihood term p(P|i3, Mi). The global 
likelihood is the normalisation term p(P|Mi), which is equal 
to the total probability: 



p{V\Mi 



p{B, ^Mi)p{V\B, Mi)d$di3. (6) 



by marginalising the joint posterior probability: 

p{B\V, Ml) = j p{B, Mi)d$. (7) 

We therefore ensure that an adequate B possesses Stokes V 
profiles at some phases that match all the observations. From 
this B posterior probability density, we can then explore the 
plausible region of the parameter space of an oblique dipole 
field given the data, assuming that the model is true. 

In the case of a limited number of spectropolarimet- 
ric observations, we will be generally interested in the field 
strength value, as only weak constraints can be placed on 
the geometrical parameters. This is particularly relevant in 
the case of a non-detection, where we wish to put an upper 
limit on the strength of an undetected dipole field. To ob- 
tain the posterior probability density for a given parameter, 
we need to marginalise the joint probability p{B\T>^ Mi) over 
the other parameters. For example, the posterior probability 
density marginalised for the field strength is: 



/ 



p{Bp\V,Mi)= p{B\V,Mi)didl3. 



(8) 



Another powerful application of Bayesian statistics is 
the ability to quantitatively test the plausibility of one hy- 
pothesis versus another. We can therefore compare the plau- 
sibility of the oblique dipole model Mi, by computing the 
so-called odds ratio of this model compared to the model 
Mo representing the absence of a magnetic field. To do so, 
we need to compute the posterior probability p{Mi\D,I) of 
a given model which is, according to the Bayes theorem: 

p{M,\I)p{D\M,,I) 



p{M,\DJ): 



(9) 



p{D\I) 

The prior term p{Mi\I) encodes the plausibility of the 
model, given our current knowledge, and the global likeli- 
hood p{D\Mi, I) encodes the plausibility of the model, given 
the new observations. 

Therefore, the odds ratio of our two competing models. 
Mo and Mi, can be written as: 

p{Mo\V,I) _p{Mo\I)p{V\Mo,I) 



odds (Mo /Ml 



(10) 



We then treat the set of phases $ as nuisance parameters 



p(Mi|P,/) p(Mi|/)p(P|Mi,/)- 

Typically, as no model is preferred prior to the acquisition of 
Stokes V observations, the ratio of priors p{Mo\I) /p{Mi\I) 
will be set to unity. The global likelihood of the dipole model 
is given by Eq. (|6|. As the model representing the absence 
of magnetic field has no parameters, its global likelihood is 
simply the product of the likelihoods that V = for each 
observation. According to [Jeffreys] p998), the evidence in 
favour of a model is considered weak when the odds are 

> 10°-^ (^^3:1), moderate when > 10^ (-^10:1), strong when 

> 10^-^ (-30:1) and very strong when > 10^ (-100:1). 



3.1 Practical implementation 

Given the complexity of disk integration for a rotating mag- 
netic star, it is not practical to solve the Bayesian problem 
analytically (for an example applied to the local Stokes pro- 
files of the Sun, see Asensio Ramos 2 011| ). A Markov-chain 
Monte-Carlo method works by choosing series of parame- 
ter sets, the parameter regions with a high likelihood being 
more likely to be picked. The sampling itself therefore re- 
flects the posterior probability density. This method is not 
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well suited to this problem because a large number of calcu- 
lations are required to assess the probability of each B config- 
uration, due to the marginalisation over all possible phases. 
We therefore chose a numerical, brute-force approach in or- 
der to explore the parameter space by sampling it with a 
regular grid. 

When evaluating Eqs. ^ and ([7|, we can make 
use of the fact that the prior probability densities for 
each parameter are independent such that $|Mi) = 
p(i3|Mi)p($|Mi). Furthermore, the prior probability for 
the phases are the same such that p{(pn) = Fi- 
nally, the likelihood $,Mi) can be factored into 
Iin=iP{Dn\B,(pnj Ml) givcu that the likelihood of one ob- 
servation Dn is only dependent on the phase cpn- Therefore, 
Eq. |7|can be written as: 



p{B\D,Mi)=p{B\Mi) 



U^^Jp{mi)p(Dn\B,^,Mi)d^ 
p{V\Mi) 



(11) 

For each point in a [Bp, i, /3, cp] parameter space, we can 
compute the likelihood between the synthetic profile pro- 
duced by these parameters and each observation. In fact, 
given the Stokes V profile symmetry properties described in 
Eq. (|3|, only the synthetic profiles for a quarter of the pa- 
rameter space need to be computed. Furthermore, synthetic 
profile interpolation is possible in the Bp direction, as the 
local Stokes V profile amplitude increases linearly although 
the shape remains the same. The maginalisation integrals 
are performed numerically using a five point Newton- Cotes 
algorithm. An adequate sampling of the Mi parameter space 
was chosen to obtain accurate values of p(Mi|P,/). 



3.2 Probability densities 

For each parameter of the oblique dipole model, a prior 
probability density is needed to evaluate Eq. ^ . The dipole 
field strength Bp is a parameter that can vary over several 
decades (from a few dozen G to a few kG). We therefore 
used a Jeffreys prior, which sets an equal probability per 
decade and is therefore scale invariant. We used a modified 
form ( Gregory||2005a ) in order to eliminate the singularity 
at B^ = 0: 



p{Bp\Mi) = 



1 



(12) 



This prior function behaves as a fiat prior when Bp < a, and 
as a Jeffreys prior when Bp > a. The maximum dipolar field 
strength Bp, max is adjusted according to the strength of the 
Stokes V signal and by the quality of the data (for example, 
if the s/n is lower, a field of a higher strength can hide in 
the noise in the case of a non-detection). For this paper, the 
grid has 250 Bp values. A finer sampling might be required 
when using a large number of observations, in order to avoid 
under sampling the joint posterior probability density. We 
tested the effect of the choice of the a parameter, and we 
settled using a value about twice the Bp grid step. 

The inclination of the rotation axis to the line of sight 
is a "position" parameter (invariant under a shift of zero po- 
sition), for which a fiat prior would be well suited. However, 
we know that if the orientations of the rotation axes of stars 



are generally randomly distributed, the probability that the 
inclination angle is between i and i -\- di is: 



p{i\Mi) = ^ sini. 



(13) 



The inclination could in principle vary from 0° to 180°. How- 
ever, if a non-null vsini is measured, it is possible to put a 
lower limit on the inclination by estimating the break-up ve- 
locity - typically around 700kms~^ for a massive OB star. 
We used a grid of 37 i values, giving a sampling at least 
every ^5°. 

The same reasoning could apply to the obliquity angle 
/3 of the magnetic axis to the rotation axis. However, we 
ignore at this point if the magnetic axes are generally ran- 
domly distributed, or if a general relation exists for the po- 
sition of the magnetic axis relative to the rotation axis. For 



example, Landstreet & Mathys ( 2000 ) have shown that the 



magnetic axis of long-period ApBp stars is generally aligned 
with their rotation axis. Therefore, instead of using a prior 
that assumes a randomness of the magnetic axis position, 
we kept a fiat prior for the obliquity angle: 



pWIMi 



1 



/3n 



(14) 



The obliquity varies from /3min = 0° to /3max = 180°, sam- 
pled every 5°. 

A fiat prior probability is also used for the rotational 
phase: 



p{ip\Mi) = 



(15) 



"^max "^min 

with (fmin = 0° and (fma^ = 360°, sampled every 5°. 

The likelihood function of a given data set Dn corrupted 
with Gaussian noise is given by: 



p{Dn\B,ip,Mi) = (27r)-^ 



exp 



(16) 



where dk represents one of the M datapoints with its vari- 
ance cTfc, and fk represents the model prediction. 

When performing a parameter estimation, assuming the 
veracity of a model, it is a good practice to treat the variance 
of our data as a model parameter itself, using our estimate 
of the variance (i.e. the error bars) as a starting point. Pro- 
ceeding in this way, the final results will be less sensitive to 
our estimation of the variance, as noise propagation can be 
problematic in heavily processed data, as can be the case 
for LSD profiles. Furthermore, this approach treats any fea- 
tures that cannot be explained by this particular model (for 
example higher order polar field components or distortions 
of the profile due to abundance spots) as additional noise. 
According to the maximum entropy principle, a Gaussian 
distribution will be the most noncommittal about informa- 
tion we do not have, leading to the most conservative esti- 
mates of the parameters. If our estimate of the variance of 
each point is denoted by Sfc, we introduce a "noise scaling" 
parameter h ( [Gregory|2005b| ) , in order to estimate the total 
noise variance Gk'- 



2 



^2 



(17) 
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where b varies around unity. The resulting expression for the 
hkehhood is: 



p{Dn\B,ip,b,Mi) = (27r) 



exp 



(18) 



The resulting posterior probability density will therefore be 
marginalised over b. This procedure is only used for the pa- 
rameter estimation of model Mi. The noise scaling param- 
eter is a scale parameter, and a Jeffreys prior will be used: 



p(6|Mi 



bin 



(19) 



The noise scaling parameter for the parameter estimation 
varies from 6min = 0.1 to 6max = 2, corresponding to an 
underestimation of the variance by a factor of 10 and an 
overestimation of the variance by a factor of 2, respectively. 
This conservative range can be verified a posteriori, by the 
probability density marginalised for b. 



4 NUMERICAL TESTS 

In order to demonstrate our Bayesian method, we have sim- 
ulated some realistic synthetic observations, and analysed 
them using our procedure. 



4.1 Ideal non detection 

We start with a simple test representing an ideal non- 
detection. The local intensity line profiles are represented 
by Gaussians, with the width of the point spread func- 
tion of an instrument with spectral resolution of 4.2kms~^ 
(R=65000), typical of current high-resolution spectropo- 
larimeters like ESPaDOnS, and an extra turbulent broad- 
ening of 5kms~^, typical of the broadening encountered in 
hot stars. We used a shallow line with a depth of 0.17 /c. The 
rotational velocity field is set to v sin i 54kms"^ The as- 
sociated circular polarisation is set to zero (i.e. no magnetic 
field) and the error bars are set to represent a certain s/n. 
We used two simulated observations corresponding to a s/n 
of 12 000 and 24 000 respectively, values that can typically 
be achieved for the LSD profiles of hot stars. The effective 
Lande factor and the wavelength of the spectral line are set 
to 1.2 and 5 000 A respectively. We performed the Bayesian 
analysis on a grid extending up to IkG, with a maximum 
rotation velocity of 700kms~^, leading to a minimum incli- 
nation of 8°. 

Figure |4] shows the probability densities resulting from 
our Bayesian analysis, for the s/n of 12 000. The three 
bottom panels show the posterior probability densities 
marginalised for each parameter - Bp, i and /3, from left 
to right respectively. The densities have been normalised 
by their maximum, in order to facilitate the display. The 
three top panels show the 2D posterior probability densities, 
marginalised for the Bp-/3, Bp-i and i-/3 planes. The contours 
encircle the regions that contain 68.3%, 95.4%, 99.0% and 
99.7% of the probability. 



135 
90 
45 
180 
135 
90 
45 







(CT ■ 






V , , , 







250 500 750 1000 45 90 135 45 90 135 180 



, (G) 



i (°) 



iS (-) 



Figure 4. Parameter estimation for the ideal non-detection 
(low s/n). Bottom: Posterior probability density functions (PDF) 
marginalised for the dipole field strength, the rotational axis incli- 
nation and the magnetic obliquity, from left to right respectively. 
The PDFs have been normalised by their maximum. Top: 2-D 
posterior probability density marginalised for the Bp-i, Bp-(3 and 
i-(3 planes. The 68.3%, 95.4%, 99.0% and 99.7% credible regions 
are shaded in dark to pale colours respectively. 



As it can be seen from the Bp probability density distri- 
bution, the bulk of the probability is concentrated at low Bp 
values. The amplitude of the Stokes V profiles produced for 
low Bp values will remain under the noise level, for all possi- 
ble dipole orientations. As Bp increases, the amplitude of the 
Stokes V profiles will progressively grow over the noise level, 
and only a restricted number of possible orientations will 
produce Stokes V profiles with amplitudes below the noise 
level, as explained in section |2.2| The probability density 
marginalised for Bp therefore has a decreasing exponential- 
like shape, with a tail extending far away from the mode of 
the distribution. 

The shape of the 2D Bp-i and Bp-^ planes gives us in- 
sight about the geometries that are more likely to produce a 
non-detection for a large field strength value. When we look 
at a star rotational equator-on {i = 90°), the obliquity re- 
quired to obtain a low amplitude Stokes V profile is around 
0° or 180°. Thus, as the field is nearly aligned with the rota- 
tional axis, the amplitude of the Stokes V profiles will stay 
small for the whole phase range. For lower inclination angles, 
the Stokes V signal will vanish only when the dipole field is 
seen nearly magnetic equator-on, therefore such an null ob- 
servation is only possible on a restricted interval of rotation 
phases (as it was shown in Figure |3|. Consequently, if we 
assume that a strong field is present but is not producing 
any Stokes V signal, it is more likely that we are observ- 
ing a dipole configuration for which the Stokes V profiles 
always have a small amplitude, rather than a dipole con- 
figuration for which the Stokes V signal vanishes for only 
a small fraction of the rotational period. This explains why 
the probability density is more extended toward higher Bp 
values for i = 90° and /3 = 0°, 180°. 

The 2D i-f3 plane shows that we have no constraint 
on the inclination and obliquity angles. The shape of the 
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Table 1. Base ten logarithm of the odds ratio log(Mo/Mi) for 
the ideal non-detections. 



Test 


log(Mo/Mi) 


log(Mo/Mi) 


log(Mo/Mi) 




low s/n 


high s/n 


joint obs. 


(1) 


(2) 


(3) 


(4) 


Flat prior 


0.976 


1.267 


1.355 


Used prior 


0.317 


0.444 


0.475 



probability density marginalised for i is in part due to the 
extra probability at higher Bp for i ^ 90°, but is mainly 
driven by the prior probability that favours large inclinations 
(as given in Eq. |T3|). Following the same reasoning, the 
aligned dipoles (/3 ^ 0° or 180°) are slightly favoured. 

Table [l] compiles the base ten logarithm of the odds ra- 
tio computed for each observation taken individually, and 
for the two observations combined together. In order to il- 
lustrate the effect of the priors, we also computed the odds 
ratios while using a flat prior for all the parameters. For a 
single low s/n observation, the model Mo for the absence of 
a magnetic field is preferred by a factor of 2. When doubling 
the s/n, the odds in favour of Mo only increase to a factor 
2.8. When the two observations are combined together, the 
odds ratio goes up to a factor of 3. 

If we had used a flat prior for all parameters, the odds 
ratios would have been in favour of the model Mo by a factor 
of roughly 10 and 20, for the low and high s/n respectively. 
Jeffreys priors are used when we ignore the exact scale of 
a parameter. This serves to balance the high scales and the 
lower scales. In a non-detection case, the bulk of the proba- 
bility is situated at low Bp values. With a Jeffreys prior, the 
small scales have as much weight as the larger scales, hence 
the configurations at high Bp values that produce bad fits 
dominate less the global likelihood. The flat prior effectively 
gives more weight to the larger scales, hence rejecting more 
strongly the Mi model. We note here that the Mo model 
is in fact a subset of the Mi model, given that the field 
strength range extends down to OG. This means that both 
models can reproduce this simulated dataset perfectly. How- 
ever, the Ml model is penalised by its complexity, which 
is mathematically encoded by the priors. Therefore, for a 
non-detection, the odds ratios will never be very strongly in 
favour of the Mo model by a large number, as the discrimina- 
tion comes mainly from "Occam's razor" . Furthermore, the 
posterior probability density is sensitive to the exact choice 
of the prior in this case. Given the decreasing exponential 
behaviour of the probability density, we think that the more 
conservative Jeffreys prior is more suitable, as it is less eager 
to reject the Mi model. 

Assuming that the dipole model is true, we can per- 
form a parameter estimation to determine which dipole 
strengths would be admissible by our observations. This 
can be achieved using the marginalised probability density 
for Bp. Integrating this probability density between two Bp 
values gives the probability that the true value of the field 
strength lies between these two values. The probability den- 
sity itself can be used in many applications, such as building 
a statistical field strength distribution for a stellar popula- 
tion. 

In other applications, such as the study of magnetically 
confined wind shocks, it would be valuable to get an up- 
per field strength limit. When a probability density has a 



Table 2. Credible region upper limits for the ideal non-detection. 



Test 


68.3% 


95.4% 


99.0% 


99.7% 




(G) 


(G) 


(G) 


(G) 


(1) 


(2) 


(3) 


(4) 


(5) 


low s/n 


30 


112 


258 


456 


high s/n 


18 


66 


157 


297 


Joint 


16 


51 


110 


196 


Joint, flat prior 


29 


138 


398 


734 


Joint, no scale noise 


20 


67 


142 


251 



Gaussian-like shape, it is customary to express our confi- 
dence intervals by regions enclosing a certain percentage of 
the total posterior probability, namely 68.3%, 95.4% and 
99.7%. With a Gaussian distribution, the 95.4% region will 
be twice as extended as the 63.8% region, and the 99.7% will 
be 3 times as extended as the 68.3% region, analogous to the 
1, 2, 3 cr contours in frequentist statistics. However, as our 
probability distribution does not have a Gaussian shape, but 
is rather shaped like a decreasing exponential, these credible 
regions will not be regular (i.e. the 99.7% credible region will 
reach much farther than 3 times the 68.3% credible region, 
taking into account the extended tail of the distribution). 
We therefore added a credible region enclosing 99.0% of the 
probability. 

Table [2] gives the upper limits of the credible regions 
(the lower limit being G in each case) for each percentage 
threshold of the probability. For example, the probability 
that the field strength is lower than 258 G is 99.0% for the 
low s/n case. For the higher s/n case, all the credible re- 
gions are narrower. When we combine the two observations 
together, the 68.3% and 95.4% regions are similar to the high 
s/n case, but the probability density is more peaked (i.e. it 
has a less extended tail toward the high Bp values), as can 
be seen from the narrower 99.0% and 99.7% credible regions 
(110 G and 196 G respectively). When combining two non- 
detection observations, the probability density of high Bp 
values narrows around i — 90° and /3 = 0°. The likelihood 
of the i = 90°, /3 = 0° configurations stays the same, as none 
of the phases produce any Stokes V signal. However, the like- 
lihood of the other degenerate configurations decreases, as 
it becomes less likely that we have observed the star both 
times during the same narrow range of phases that do not 
produce a Stokes V signal, if the observations were taken 
at random times. Therefore, the high- Bp tail becomes less 
important as we combine multiple observations. 

Table [2] also gives the credible regions for the combined 
observations, while considering a flat prior for each parame- 
ter. In that case, the shape of the probability density is more 
weighted toward the high scales, and the credible regions are 
more extended toward the high Bp values. 

We also give the credible regions when considering a 
fixed variance of our data. However, in this particular case, 
as the observations can be reproduced perfectly by both 
models (as we set Stokes V strictly to zero), the inclusion 
of the noise scaling parameter does not change the prob- 
ability density much, except for a slight tightening of the 
credible regions as we allowed for the possibility of variance 
overestimation. 
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4.2 Ideal detections 

In order to show the behaviour of the algorithm in the pres- 
ence of a magnetic signal, we now add non-zero Stokes V 
profiles to the simulated observation corresponding to a s/n 
of 12 000. Once again, no random noise was added, to illus- 
trate the ideal case. To draw a connection between the cred- 
ible regions found in the preceding section, we chose field 
values corresponding roughly to the upper limits of these 
credible regions: 30 G, 100 G, 250 G and 450 G. The mag- 
netic configuration was set to i = 90° and /3 = 90°, at the 
rotational phase when we are looking straight at the mag- 
netic pole (cp — 0°). This corresponds to a profile with a 
shape like the leftmost profile of Figure [2] 

Table [3] gives the odds ratios for the used priors, as 
well as the ones computed with flat priors for comparison. 
We also give the credible regions containing 68.3%, 95.4%, 
99.0% and 99.7% of the total probability. The lower limit is 
OG, unless indicated otherwise. 

The posterior probability density for the 30 G observa- 
tion is indistinguishable from the V — ^ observation from 
the previous section. The odds ratio is still marginally in 
favour of the Mo model, and the credible regions for Bp are 
similar. 

At 100 G, the odds ratio is still in favour of Mo, although 
both models are now nearly as likely (log(Mo/Mi = 0.097)). 
The odds ratio computed with a flat prior still rejects the 
Ml model by a factor of three. The shape of the probability 
density, as shown in Figure |5] is also different than the non- 
detection case, with a secondary peak in the probability den- 
sity marginalised for Bp, although the probability density 
still peaks at OG. The credible regions for Bp are therefore 
extended toward higher values, and the 68.3% credible re- 
gion upper limit (123 G) encompasses the real field value. As 
it was the case for the ideal non- detection, an aligned mag- 
netic configuration is preferred, as it maximises the chances 
of observing this particular profile shape. The 2D plane 
is different from that of the ideal non-detection. The shape 
of the magnetic signal adds an extra constraint, by reject- 
ing the configurations for which the positive pole is never 
located on the visible hemisphere. 

At 250 G, the odds ratio has switched in favour of the 
dipole model, by three orders of magnitude (log(Mo/Mi) = 
—2.98). However, our chosen prior is again more conserva- 
tive than the flat prior, which favours Mi more strongly 
(log (Mo /Ml) = —3.22). As a flat prior overweights the 
larger scales compared to a Jeffreys prior, our used prior 
needs a more significant likelihood at large Bp values in or- 
der to favour the magnetic model. 

The posterior probability density is now typical of a de- 
tectable magnetic field (Figure |6|. A sharp cut in the prob- 
ability density marginalised for Bp and the 2D planes shows 
a tight constraint on the lower field limit. The Bp distribu- 
tion peaks around 275 G and decreases slowly toward higher 
Bp values. Therefore, the credible regions for Bp are not 
symmetric with respect to the mode of the distribution. 

The lowest field strength values correspond to configu- 
rations for which the positive magnetic pole is located at the 
stellar disk center. Larger field values are possible when we 
are looking at the positive pole from an angle. For example, 
if the field is aligned (/3 = 0°), the amplitude of the Stokes 
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Figure 5. Same as Figure^ for an ideal detection with a field 
strength of 100 G. 



V profiles will decrease as the inclination decreases, until we 
reach i ^ 0° where the magnetic signal vanishes. 

Because there is only one observation, the aligned con- 
figurations are preferred, as shown by the probability density 
marginalised for /S. Obviously, when the inclination is 90°, an 
aligned configuration would not provide any Stokes V signal 
and cannot reproduce the observed profile. Some obliquity 
would be necessary in order to put the positive pole on the 
visible hemisphere, which explains the dip at z = 90° in the 
probability density marginalised for i and marginalised for 
the Bp-i plane. 

When i ^ 0° and the field is aligned, we are always look- 
ing directly at the positive pole, which makes this configura- 
tion quite likely. However, our prior knowledge tells us that 
a low inclination is less likely than a high inclination. For 
this reason, the posterior probability density marginalised 
for i decreases toward i ^ 0° and 180°. 

At the present B configuration {i = 90° and /3 = 90°) 
the particular profile shape (anti-symmetrical) only occurs 
during a certain phase range (see Figure |2| , which makes 
such a configuration less likely than the aligned configura- 
tions. Therefore, with only one observation, no meaning- 
ful constraints can be put on the angles. However, the field 
strength can be better constrained, and the probability den- 
sity marginalised for Bp can provide statistical insight into 
the more likely field strength values. 

At 450 G, the odds ratio favours Mi by more than 13 
orders of magnitude. When the Mo model is so strongly 
rejected, the difference in odds ratios between our chosen 
prior and the flat prior is less important than when the signal 
detection was more marginal. This is because the bulk of 
the likelihood is located at higher Bp values, where the two 
priors are similar. The posterior probability density is quite 
similar in shape to the one with a 250 G field, with the Bp 
distributions shifted to higher values, as it can be seen from 
the credible regions in Table |3] 
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Table 3. Odds ratios (log(Mo/Mi)) for the adopted and flat priors, as well as credible regions 
(lower limit of G unless indicated otherwise) for the ideal detections. 



Test 




log(Mo/Mi) 


log(Mo/Mi) 


68.3% 


95.4% 


99.0% 


99.7% 






Used prior 


Flat prior 


(G) 


(G) 


(G) 


(G) 


(1) 




(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


30G ( - 


68.3%) 


0.303 


0.939 


34 


128 


290 


498 


lOOG ( - 


95.4%) 


0.0970 


0.476 


123 


338 


625 


825 


250G ( - 


99.0%) 


-2.98 


-3.22 


196 - 407 


155 - 763 


143 - 938 


141 - 991 


450 G ( - 


99.7%) 


-13.4 


-13.9 


410 - 658 


384 - 928 


379 - 995 


363 - 1000 
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Figure 6. Same as Figure [4] for an ideal detection with a fleld 
strength of 250 G. 



4.3 Realistic case 

So far, we have explored the behaviour of the algorithm in 
cases where the models can reproduce the data perfectly. 
Obviously, real observations will have some deviations from 
the predicted values, introduced by the noise. We will now 
explore how the noise corruption affects our detection ca- 
pacity. 

The top spectrum of Figure [7| shows a simulated ob- 
servation, consisting of only random noise generated from a 
normal distribution corresponding to a s/n of 12 000. Figure 
|8] presents the posterior probability density for that specific 
observation. Notice how the Bp distribution looks like the 
probability density for the 100 G observation of the previous 
section (Figure [g]) In fact, the odds ratio is even in favour 
of the dipole model (log(Mo/Mi) = -0.419), even though 
the signal is pure noise. This is because the observation is 
better fitted by a non-null magnetic field. The absolute best 
fit to the data is given by the maximum of the likelihood. We 
illustrate this fit by the red profile overplotted on the data 
(Figure [7|. Therefore, the detection of a real signal embed- 
ded in the noise is ambiguous, as the noise could have - by 
chance - the shape of a magnetic profile. An observer would 



^ The (5 density probability is different from the 100 G ideal de- 
tection because the simulated noise pattern has a shape similar to 
a magnetic fleld whose pole is located at a non-null rotational ra- 
dial velocity. Such a location requires a dipole that is not aligned. 



Table 4. Odds ratios (log(Mo/Mi)) for the realistic simulated 
observations for the pure noise case and the Bp = 125 G case. 
The rotational phases, chosen randomly, are also given. 

TbSt log(Mo/Mi) ^ log(Mo/Mi) 

Noise only Bp = 125 G 



(1) 


(2) 


(3) 


(4) 


Obs. 1 


-0.419 


0.76 


-4.04 


Obs. 2 


0.306 


0.43 


-0.412 


Obs. 3 


0.295 


0.40 


-0.484 


Obs. 4 


0.178 


0.20 


-2.26 


Obs. 5 


0.281 


0.60 


-0.388 


Combination 


0.295 




-9.05 



likely not feel confident to report a magnetic detection based 
on only one observation like this one. 

In order to verify if this ambiguity can be lifted by re- 
observing the star, we generated 4 additional simulated ob- 
servations, again with pure normal noise (rest of Figure |7|. 
Each of these observations leads to an odds ratio in favour 
of Mo by about a factor of two (Table [4] column 2). 

The best fit achievable for each observation taken in- 
dividually is again overplotted in red. However, nothing re- 
stricts the magnetic configuration B to be the same for each 
observation. The best fit produced by a single B configura- 
tion is given by the maximum of the joint likelihood, illus- 
trated in green. Although the data can be reproduced by 
a non null magnetic field, the odds ratio of the combined 
observations is in favour of the Mo model (log(Mo/Mi) = 
0.295). The slight improvement of the fit produced by the 
oblique dipole model (see the reduced indicated in Figure 
|7|, is not enough to justify the use of a more complex model. 

The posterior probability density (Figure[9]) is similar in 
shape to the perfect non-detection that was shown in Figure 
|4] The high- Bp tail is less extended that in the case of a sin- 
gle observation. This can by seen more easily by comparing 
the 2D probability densities for the Bp- {3 and Bp-i planes 
of Figures [9] and |4] As explained in the previous sections, 
combining multiple observations has this effect because, in 
the case of non-detections, the high-i^p inclined dipole con- 
figurations become less likely, as this would mean that all 
the observations were taken during a narrow range of phase. 
Obviously, the possibility that the observations were taken 
at the same rotational phase can generally be assessed by 
the time span of the observations and the possible range of 
rotational periods. In the case where it is known that the 
rotational period is quite long compared to the observation 
time span, more information is available than is assumed 
in this algorithm (that the phase of each observation can 
assume any value). In that case, it would be wiser to com- 
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Figure 7. Realistic simulated observations of pure Gaussian 
noise. The intensity line profile is shown at the bottom, and the 
five noisy Stokes V profiles are shown at the top. The dotted 
lines display the range of the fits. The no magnetic field model 
Mq (V = 0) is shown by the black lines. The best fits of the 
oblique dipole model Mi for each observation taken individually 
(represented by the maximum of the individual likelihoods) are 
shown in red. The best fit for all the observations taken together 
by a single B geometry (the maximum of the joint likelihood) is 
shown in green. The corresponding reduced x^s are given on the 
right side, with matching colours. 
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Figure 8. Same as Figure^ for Obs. 1 of the realistic simulated 
observations of pure noise. 
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Figure 9. Same as Figure^ for the combined realistic simulated 
observations of pure noise. 



bine all the observations together and treat them as a single 
observation in order to get a more meaningful probability 
density. 

Table [s] compiles upper limits to the 68.3%, 95.4%, 
99.0% and 99.7% credible regions, extracted from the pos- 
terior probability density marginalised for Bp. Although we 
combined five observations, the 68.3% credible region for Bp 
extends to 38 G, which is higher than the value obtained for 
one perfect observation of low s/n (30 G), because of the 
deviations introduced by the noise. However, as mentioned 
above, the Bp distribution of the combined observations does 
not extend as far than as that from a single perfect observa- 
tion, as illustrated by the 99.7% credible region that extends 
only to 370 G, compared to 456 G for the perfect observation. 
We also compiled the credible regions we would obtain with- 
out the use of the noise scaling parameter. As expected, the 
deviations from the model are consistent with our estima- 
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Table 5. Credible region upper limits for the combined realistic 
simulated observations of pure noise. 



Test 


68.3% 


95.4% 


99.0% 


99.7% 




(G) 


(G) 


(G) 


(G) 


(1) 


(2) 


(3) 


(4) 


(5) 


With noise scaling 


38 


108 


214 


370 


Without noise scaling 


38 


107 


213 


368 



tion of the variance (the error bars) and the credible regions 
are nearly the same in both cases. 

We have therefore demonstrated that a suspicious sig- 
nal that has a shape similar to that of a magnetic signal 
can be verified by the acquisition of a small number of ad- 
ditional observations. We now demonstrate that the same 
is true for a real signal that is sufficiently embedded in the 
noise to render the odds ratio of a single observation am- 
biguous. Given that in the ideal detection case we were able 
to detect a field with a strength near the 95.4% upper limit 
of the ideal non-detection, we chose a field strength close to 
the upper limit of the 95.4% credible region of the previ- 
ous example. The B configuration is given by Bp — 125 G, 
i — 90° and /3 = 90°. Five phases were randomly gener- 
ated {lp = 0.76, 0.43, 0.40, 0.20, 0.60), and the correspond- 
ing Stokes V profiles were added to the previous simulated 
dataset of pure noise. The resulting simulated observations 
are shown in Figure [To] The underlying real Stokes V profile 
is shown as the black dashed curve for each observation. 

Table[4] gives the odd ratio for each observation (column 
4), all of which favour the dipole model Mi, although obser- 
vations 2, 3 and 5 less strongly than observations 1 and 4. In 
fact, the odds ratios of the formers are similar to the odds ra- 
tio of the first observation of the pure noise test and on their 
own, each of these observations would result in an ambiguous 
signal detection. However, combining all the observations to- 
gether, the odds ratio is now strongly in favour of the dipole 
model by 9 orders of magnitude (log (Mo /Mi) = —9.05). 

Figure [TT] shows the posterior probability densities for 
the combined observations. There is a sharp lower limit on 
the possible magnetic strengths, illustrated by the probabil- 
ity density marginalised for Bp, and for the 2D B^-i and 
Bp-P planes. Given that most of the likelihood is situated 
at non-null field strengths, the improvement of the fit to the 
data justifies the m ore c omplex model, as shown by the re- 
duced on Figure 



10 



for the Mo model fits (black lines) 
and the fit produced by maximum of the joint likelihood for 
Ml (green curves). 

Different choices are possible in order to express the 
derived values for the model parameters. One can choose 
for example the maximum of the joint posterior probability 
density (MAP), or the mode of the marginalised probabil- 
ity density for each parameter. Note however that if the 
probability distribution is complex, the parameters given by 
the MAP do not necessarily correspond to the mode of the 
marginalised probability densities. Usually, the MAP will 
produce the best fit to the data given that the a-priori in- 
formation does not exclude any interesting parts of the pa- 
rameter space, but does not necessarily represent the bulk 
of the probability. The mode of each parameter represents 
well the bulk of the probability, but does not necessarily 
give an excellent fit to the data. Using the median of the 
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Figure 10. Realistic simulated observations of noise with a 
Stokes V signal corresponding to Bp = l2hG (i = (3 = 90°). The 
intensity line profile is shown at the bottom, and the five noisy 
Stokes V profiles are shown at the top. The dotted lines display 
the range of the fits. The underlying real Stokes V signals are 
shown with black dashed curves. The null magnetic field model 
{V = 0) is shown by the black lines. The best fit by a single B 
geometry for all the observations taken together (the maximum 
of the joint likelihood) is shown in green. The fit produced by the 
MAP B parameters is shown in blue, and that produced by the 
modes of the posterior probability density marginalised for each 
parameter is shown in magenta. The corresponding reduced x^s 
are given on the right side, with matching colours. 
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Table 6. Estimation of the B parameters for the reahstic simulated observations of noise plus 
a signal corresponding to Bp = 125 G. The first three columns give the maximum of the joint 
posteriori probability density as well as the modes and the medians of the joint posteriori probability 
marginalised for each parameter. The last four columns give the credible regions of each parameter. 



Parameter 


MAP 


Mode 


Median 


68.3% 


95.4% 


99.0% 


99.7% 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


Bp (G) 


210 


195 


220 


155 - 261 


118 - 508 


109 - 793 


102 - 892 


in 


61 


85 


71 


62 - 118 


38 - 142 


27 - 152 


21 - 158 




119 


95 


109 












70 


90 


66 


54 - 125 


23 - 157 


14 - 165 


12 - 168 




110 




113 











135 
90 
45 
180 
135 
90 
45 

1.8 

0.8 
0.6 
0.4 
0.2 
0.0 






Table 7. Traditional field diagnostics applied to our numer- 
ical tests. The global longitudinal field was integrated from 
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Figure 11. Same as Figure^for the combined realistic simulated 
observations of noise plus a signal corresponding to Bp = 125 G. 



marginalised probability densities is usually a good compro- 
mise (Gregory 2005b). 

In Table [6] we list the MAP, the mode and the median 
for each B parameter, which in this case are quite similar. 
The MAP, the modes and the medians all yield similarly 
good fits to the data (the MAP and the modes fits are shown 



in Figure 10 in blue and magenta respectively). 

The range of the credible regions of the probability den- 
sity marginalised for each parameter are compiled in Table|6] 
Using the median with the 68.3% credible region, we would 
infer a dipole strength Bp — 220 -tl G, which is slightly 
higher than the input dipole field strength of 125 G. We 
verified that this difference was due to the particular noise 
pattern used in the data, as additional noise simulations al- 
lowed us to recover the input value within the 68.3% region. 
The real values of the i and /3 angles are recovered by the 
68.3% credible regions, although the constraints are poor, 
as expected. 



4.4 Comparison with traditional diagnostics 

With low-resolution instruments (such as FORSl and 
F0RS2 at the Very Large Telescope), one is generally only 
sensitive to the global longitudinal component of the mag- 
netic field, as the rotationally-broadened spectral lines are 
not resolved. This global longitudinal field value is extracted 
from the spectrum using the "slope" method as described in 



1 to +55kms-^ 


Test 


Bi 


\Bi/a\ 


P{V) 




(G) 




(%) 


(1) 


(2) 


(3) 


(4) 




Ideal detections 




30 G 


11±33 


0.34 





100 G 


37 ±33 


1.12 





250 G 


93 ±33 


2.76 


52.0 


450 G 


168 ±33 


5.03 


99.99997 




Noise only 




Obs. 1 


72 ±33 


2.16 


60.8 


Obs. 2 


12 ±33 


0.38 


64.1 


Obs. 3 


2 ±33 


0.08 


16.7 


Obs. 4 


-16 ±33 


0.48 


87.9 


Obs. 5 


10 ±33 


0.32 


13.6 




Bp = 


125 G 




Obs. 1 


74 ±33 


2.22 


99.5 


Obs. 2 


-28 ±33 


0.87 


91.9 


Obs. 3 


-35 ±33 


1.07 


64.8 


Obs. 4 


0±33 


0.01 


99.8 


Obs. 5 


-27 ±33 


0.82 


54.7 



Bagnulo et al. (2002 2006 ). The presence of a magnetic field 
in a single observation is therefore diagnosed by the signifi- 
cance of the global longitudinal field measurement compared 
to its error bar. 

In the case of high-resolution spectropolarimetry (e.g. 
ESPaDOnS at Canada- France- Hawaii Telescope, Narval at 
Telescope Bernard-Lyot or HARPSPol at ESQ La Silla 3.6 m 
Telescope), the rotationally-broadened spectral lines are re- 
solved and the field is generally diagnosed by the deviation 
of the circular polarisation profile with respect to V = 0. 
The detection can be quantified by the probability that such 
a deviation from 1/ = is produced by random noise (the 
false alarm probability or FAP). The detection probability 
P can then be expressed as P = 1 — FAP. Following |Do-| 



nati et al. (1997), a field is generally considered detected if 
the FAP is less than 10~^ (P > 99.999%), and marginally 
detected when FAP is less than 10"^ (P > 99.9%). With 
high-resolution spectropolarimetry, it is also possible to in- 
tegrate the signal over the line profile in order to recover 
a value equivalent to the global longitudinal field obtained 
from low-resolution instruments ( Donati et al.|p^97 Wade 



et al.||2QQ Qa). Although it is possible to detect a magnetic 
signal even when the global longitudinal field is null, this is 
a useful quantity for producing longitudinal field curves and 
for comparing with low-resolution data. 

We applied these traditional diagnostics to our sim- 
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ulated datasets and we report the results in Table jT] In 
the case of the ideal detections, both the longitudinal field 
(columns 2 and 3) and the detection probability (column 
4) yield a detection when the field strength reaches 450 G. 
However, the odds ratios are already in favour of the mag- 
netic model at 250 G. The detection probability only looks 
at the total deviation, whereas the odds ratio looks for a 
shape similar to that of the theoretical model. 

For the realistic case consisting only of noise, all five 
observations have longitudinal field measurements below 3cr 
significance (see column 3). Moreover, no signal is detected 
in the Stokes V profiles, as shown by the detection probabil- 
ities. The same is also true for the realistic case consisting 
of noise plus a signal for Bp — 125 G, although the detection 
probability nearly reaches a marginal detection for Obs. 1 
and 4. For these observations, the odds ratios were in favour 
of the magnetic model by 2 and 4 orders of magnitude re- 
spectively. For the remaining observations, the odds ratios 
were also in favour of the magnetic model, although at a 
lower significance. Therefore, if we had observed any of the 
simulated observations with Bp — 125 G, the traditional di- 
agnostics would not have diagnosed the presence of a field, 
but the odds ratio would have indicated the possible mag- 
netic signal. A few additional observations are enough to 
distinguish between real noise and a buried signal consistent 
with an oblique dipole field, as demonstrated by the example 
here. Therefore, Bayesian odds ratios provide a quantitative 
indication of stars worth re-observing in magnetic surveys. 



5 APPLICATION TO REAL DATA 

In order to present the application of this method to real 
data, we will use high-resolution spectropolarimetric ob- 
servations of the magnetic B-type star LP Ori (=Parl772, 
HD 36982). 

LP Ori is often considered a chemically peculiar star 
because it was first classified as a B1.5p star by |Sharpless| 
(|1952V LP Ori's status as a He-strong or He- weak star is still 
uncertain. An inspection of our seven spectra (described be- 
low) does not reveal any significant variation in the He-line 
strength that would indicate a He- strong or He- weak star. 
Furthermore, comparing our spectra with the BSTAR grid 
of synthetic spectra from non-LTE tlusty models ( Lanz fc| 
Hubeny||2007J , we find a reasonable agreement with a tem- 
perature of 20 kK, a log^ of 4.0 and a vsini of 80kms~^. 
LP Ori is also a candidate Herbig Ae/Be star, because of its 
far-infrared excess, although no emission is present in the 
visible spectrum. [Manoj et al.| ( |2002| suggested that LP Ori 
is an object transiting from the pre-main sequence to the 
main sequence. LP Ori seems to be a single star. No radial 
velocity variation was found by Abt et al. ( 1991 ) nor was any 
speckle companion with a K-band ratio of less than 0.04 for 



a separation of 150 mas or more ( Preibisch et al.|1999|. 

The magnetic field of LP Ori was first reported by [Petit] 
et al.| p008). In that paper, they used three Stokes V ob- 
servations obtained with the ESPaDOnS and Narval spec- 
tropolarimeters to detect the magnetic field. 

ESPaDOnS and Narval are twin high-resolution spec- 
tropolarimetric instruments located at Canada-France- 
Hawaii Telescope and Telescope Bernard-Lyot respectively. 
A polarisation measurement consists of a set of 4 sub- 



exposures taken with different polarimeter configurations. 
From this measurement set, the circular polarisation Stokes 
V spectrum is extracted, as well as a diagnostic null polarisa- 
tion spectrum (labeled N) by combining the sub-exposures 
in such a way that the astronomical object's polarisation 
should cancel out. ESPaDOnS frames were processed us- 
ing the Upena pipehne provided by CFHT. The Narval 
frames were processed by the TBL archive pipeline. Both 



pipelines use the reduction package Libre-ESPRIT (Donati 
et al.||1997| ). The spectral range of both instruments covers 
the 370 nm to 1050 nm wavelength band, with a resolution 
i?- 65 000. 



Petit et al. ( 2008 ) used the Bayesian method described 



here to estimate the dipole strength of LP Ori, obtaining 
1150 ^200 G. For the present analysis, we have obtained one 
additional ESPaDOnS observation (within the context of the 
MiMeS CFHT Large Program) and an additional 3 Narval 
observations. The observation log, which includes the new 
observations as well as those analysed by Petit et al. (2008), 
is given in Table 



The sample of observations has a large 
range of s/n, and is therefore suitable for testing the be- 
haviour of the Bayesian method on real data. 

As is customary, we applied the LSD procedure to 
our observations. We used the iLSD code described by 
Kochukhov et al. ( 2010| ). We chose spectral lines from a Vi- 
enna Atomic Line Database (VALD; Kupka et al.|20 00) list 
corresponding to an atmosphere model with T = 20 kK and 
\ogg = 4.0. From that list, we chose metallic lines and weak 
He lines that were unblended with strong Balmer lines and 
uncontaminated by telluric lines or strong nebular emission. 
The line depths have been slightly adjusted to match the 
spectra. The final line list is shown in Table |11| The in- 
tensity (d) and polarisation (dgX) weights of each line were 
normalised by 0.2 and 120 respectively. Therefore, the dis- 
played 2/-axis scale of the resulting Stokes V and diagnostic 
null LSD profiles (Figures 12 and 13) corresponds to a line 



with a d = 0.2 and 600. 

The two sharp lines in the blue continuum of the Stokes 
/ LSD profiles are residuals from telluric lines adjacent to 
the HeIA7281 line. The emission bump present in the LSD 
profiles of the two first Narval observations are residuals 
from He emission, most likely of nebular origin because their 
centroid velocities and scaling match the nebular Balmer line 
emission. 

We computed the traditional diagnostics - the detection 
probability of the individual profiles and the global longi- 
tudinal field - by integrating zbll5kms~^ around the line 
centre, for both the Stokes V and N profiles. The results 
are displayed in Table [s] Definite detections (column 8) are 
achieved for the three ESPaDOnS observations. Marginal 
detection is achieved only for one Narval observation. No 
signal is detected in the null profiles (column 11). 



^ During the analysis of the complete set of observations of 
LP Ori for this paper, it was discovered that the sign of the Stokes 
V profiles corresponding to two spectra employed by Petit et al. 
(2008) was inverted (the Jan. 2006 ESPaDOnS spectrum and the 
Narval spectrum). During the present analysis we have corrected 
the sign of these spectra and verified the sign of all others. We 
have verified that the results of Petit et al. (2008) are not sub- 
stantially modified by this change. 
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Table 8. Log of observations for LP Ori. The exposure time is given as the total of the sub-exposures. The signal-to-noise ratio is the 
mean, per 1.8 km s"-*^ spectral pixel, between 500 nm and 600 nm. The last six columns give the global longitudinal field, the significance 
of the longitudinal field measurement and the signal detection probability for both Stokes V and the null N profiles. 



Date 


Instr. 


HJD 


texp 


s/n 


Bi V 


\Bi/a\ V 


P V 


Bi N 


\Bi/a\ N 


P N 






(+2 450 000) 


(s) 




(G) 




(%) 


(G) 




(%) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


(7) 


(8) 


(9) 


(9) 


(11) 


2007-11-08 


Nar. 


4413.54984 


4000 


73 


300 ± 434 


0.69 


60.0 


441 ±436 


1.01 


31.9 


2007-11-09 


Nar. 


4414.54299 


6000 


170 


93 ± 174 


0.53 


99.8 


-293 ±173 


1.69 


66.2 


2007-11-10 


Nar. 


4415.50631 


8800 


232 


199 ±120 


1.66 


87.1 


-3 ±120 


0.03 


96.7 


2007-11-11 


Nar. 


4416.54958 


6000 


244 


354 ±113 


3.13 


99.998 


-18±112 


0.17 


9.4 


2006-01-11 


ESP. 


3747.79588 


9600 


312 


314 ±91 


3.45 


100 


-244 ±91 


2.67 


64.0 


2007-03-06 


ESP. 


4166.76514 


9600 


510 


98 ±50 


1.93 


100 


-1±51 


0.02 


66.4 


2010-02-23 


ESP. 


5250.77084 


6400 


426 


347 ±61 


5.62 


100 


97 ±62 


1.57 


97.7 



The global longitudinal fields for the Stokes V LSD pro- 
files (column 6) show a positive trend, hinting that the posi- 
tive magnetic pole is located somewhere on the visible hemi- 
sphere. Although the longitudinal field approaches zero, the 
magnetic signal is still detectable given a sufficient s/n (for 
example 2007-03-06). This hints that we are looking nearly 
at the magnetic equator at some phases. 

Assuming a reasonable range of possible rotation pe- 
riods, defined by the vsini and breakup velocity, the lon- 
gitudinal field curve and the LSD profiles can be phased 
with various periods. Figure shows the longitudinal field 
measurements (black dots) , along with sinusoidal curves for 
three possible periods. As no other variability has been ob- 
served by other means (spectral or photometric), this is a 
good example of a case where the rotational phases of the 
observations are not known. 

Given the range of global longitudinal field measure- 
ments, we extended the grid up to Bp = 3kG. We set the 
a parameter of the Jeffreys prior to 25 G, corresponding to 
two times the parameter grid step. We performed the anal- 
ysis on both the Stokes V and null profiles. The fit of the 
Stokes / profiles is shown in Figure [12] Vertical dotted lines 
represent the velocity range of the fits. 

The odds ratios for each individual observation are 
given in Table[9] The odds ratios for the null profiles (column 
3) all favour the absence of a magnetic field (Mo model). 
The odds ratios for the Stokes V observations all favour the 
oblique dipole model (Mi), but vary from more than 10 or- 
ders of magnitude for the high s/n observations to less than 
one order of magnitude for the low s/n observations. As men- 
tioned for realistic numerical tests (Section 4.3 ), it is possible 
to obtain an odds ratio in favour of the dipole model when 
the noise happens to have a shape similar to a magnetic sig- 
nal. When combining all the observations together, we get 
an odds ratio strongly in favour of the magnetic model, by 
73 orders of magnitude, as expected given the definite sig- 
nal detections in the ESPaDOnS observations. However, if 
only the low s/n observations were available (2007-11-08, - 
09 and -10), the situation would be more ambiguous. None 
of these observations lead to a traditional definite detection 
when considered on their own. We therefore combined and 
analysed these three observations in pairs, and then all to- 
gether. The three possible pairs of observations led to odds 
ratios in favour of Mi by more than one order of magni- 
tude, and the combination of the three observations lead to 
an odds ratio log(Mo/Mi) = —3.3, i.e. more than 3 orders 
of magnitudes. Therefore, the Bayesian algorithm is able to 



Table 9. Stokes V and null profiles odds ratio for LP Ori ob- 
servations taken individually, combined together, and for various 
combinations of the low s/n observations. 



Date 


log(Mo/Mi) 


log(Mo/Mi) 




V 


N 


(1) 


(2) 


(3) 


Nar. 2007-11-08 


-0.32 


0.15 


Nar. 2007-11-09 


-1.7 


0.18 


Nar. 2007-11-10 


-0.86 


0.15 


Nar. 2007-11-11 


-11.5 


0.32 


ESP. 2006-01-11 


-13.5 


0.26 


ESP. 2007-03-06 


-24.8 


0.42 


ESP. 2010-02-23 


-16.9 


0.43 


Combination 


-73.3 


0.52 


2007-11 08±09 


-2.3 


0.20 


2007-11 08±10 


-1.3 


0.17 


2007-11 09±10 


-2.7 


0.16 


2007-11 08±09±10 


-3.3 


0.17 



recover the magnetic signal that would have been buried in 
the noise and undetected by traditional diagnostics. 

In Figures ^] and |13| the grey dashed lines represent 
the Mo model (V = 0). The best fits achievable by a dipole 
model for each observation taken individually (maximum 
of the individual likelihoods) are shown in red. The best fit 
produced by a single B oblique dipole (maximum of the joint 
likelihood) is shown in green. The reduced x^s are indicated 
with corresponding colours. For the Stokes V observations 
where odds ratios strongly favour the Mi model, the reduced 
is much improved with the addition of the more complex 
model compared to the of Mo model. Not only are 
the data reproduced by the dipolar profiles, they can all be 
simultaneously reproduced by a single B configuration, as 
shown by the similar for the individual fit (red) and the 
maximum of the joint likelihood (green). 

However, the reduced remains high in one case (2010- 
02-23; x?ed=l-95), meaning that there is some extra variance 
in the observation that neither models are able to reproduce 
(a 3(7 deviation would correspond to a reduced of 1.65). 
The reduced x^s are low for all the null profile observations 
(Figure 13). This points toward a systematic deviation or a 



model-based effect for the Stokes V observation of Feb. 2010 
rather than an extra instrumental scatter or underestimated 
error bars. 

When performing parameter estimation, the extra scat- 
ter is addressed by the noise scaling term (Eq. 17). It is 



therefore possible to extract the probability density function 
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Figure 12. Bottom: LSD line profiles of LP Ori (black). The 
synthetic line profile used for the Bayesian analysis is shown in 
red. The dashed vertical lines display the range of the fits. Top: 
LSD Stokes V profiles of LP Ori. The grey curves correspond to 
the Mo model (V = 0). The red curves correspond to the best fit 
for each observation taken individually (represented by the max- 
imum of the individual likelihood). The green curves represent 
the best fit for all the observations taken together by a single B 
geometry (the maximum of the joint likelihood). The blue curves 
correspond to the maximum of the posterior probability density 
(MAP). The corresponding reduced x^s are given on the right 
side, with corresponding colours. 



Figure 13. Same as Figure^] for LP Ori null profiles. 

marginalised for the noise scaling parameter, and determine 
if the model is able to reproduce the observations down to 
the noise level. Figure [15] shows that the Bayesian estimate 
of the variance is larger than the assumed variance (i.e. er- 
ror bars) for both Stokes V and the null profiles, but by less 
than a factor of two. 

The posterior probability density for the Stokes V ob- 
servations is shown in Figure [16] Dipole configurations with 
large inclination or obliquity are less favoured, as shown by 
the probability densities marginalised for i and /3, because 
the observations mainly show either the positive pole or the 
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Figure 14. Longitudinal field curve for LP Ori (black dots) to 
which we have superposed sinusoidal curves for three of the pos- 
sible periods. 
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Figure 16. Same as Figure^ for the combined Stokes V obser- 
vations of LP Ori. 
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Figure 15. Posterior probability density marginalised for the 
noise scaling parameter for the Stokes V (black) and null pro- 
files (dotted red) of LP Ori. 



magnetic equator on the visible hemisphere. This is more 
likely to occur if the i and /3 angles are small and the negative 
pole spends little or no time on the visible hemisphere. The 
angle values are interrelated as shown by the 1-/3 2D plane. 
If the inclination is small, the obliquity is more likely to be 
large, in order to display an equator-like magnetic signature. 
The probability density marginalised for the field strength 
shows a sharp lower limit. The high- Bp tail of the distri- 
bution is attributable to the high-inclination (low obliquity) 
configurations, as shown in the 2D planes for Bp-i and Bp-f3. 

Figure [T7| shows the probability densities for the null 
profile observations, which are, as expected, similar to the 
realistic simulation of pure noise (as shown in Figure |9| . 
We also show in Figure [18] the probability densities we ob- 
tained when considering only the three low s/n observations. 
The constraints on the parameters, especially the angles, are 
worse than when we consider the full data set. Therefore, the 
constraints on the angles are mainly defined by the high s/n 
observations. 
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Figure 17. Same as Figure^ for the combined null observa- 
tions of LP Ori. 



The credible regions of the probability density 
marginalised for Bp are given in Table [lO] We also present 
the maximum of the joint posterior probability density 
(MAP), the mode of the probability density marginalised for 
Bp, as well as the median. For the Stokes V observations, 
the MAP, mode and median values are all similar, there- 
fore the fits to the data obtained for these values are nearly 
undistinguishable. The fit obtained with the MAP values is 

Using the 68.3% credible 



shown in Figure 



12 



(blue curve) 

region and the median, the dipole field strength of LP Ori is 
estimated to be 911 ^244 G. As was seen for the probability 
density in Figure the lower limit on the field strength is 
sharp. The 68.3% credible region lower limit is 667 G, and 
557 G for the 99.7% credible region. The distribution has a 
high-Bp tail, and the 99.7% credible region extends up to 
2.6kG. 

The second row of Table [lO] illustrates the effect of the 
noise scaling parameter on the inferred field values. When 
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Table 10. Field strength estimation for Stokes V and null profile observations of LP Ori. The first three columns give the maximum of 
the joint posterior probability density as well as the mode and the median of the joint posterior probability density marginalised for Bp. 
The last four columns give the credible regions. 
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Figure 18. Same as Figure [4] for the combined low s/n Stokes V 
observations of LP Ori (2007-11-08, 2007-08-09 and 2007-11-10). 



considering the variance as a model parameter, the MAP, 
mode and median are shifted to slightly lower Bp values. 
The credible regions are somewhat larger, and also slightly 
shifted to lower values. 

The parameter estimation for the low s/n observations 
is less robust than that of the full dataset. The estimated 
field value from the median and the 68.3% credible region is 
633 ^247 G. However, we do not have a good constraint on 
the lower limit, as the credible regions quickly go to zero. 
The 99.7% credible region extends to 2.7 kG as well. 

For the null profile observations, the MAPs and modes 
are all G. Given the shape of the probability distribution, 
the median is non-null. However, given the probability den- 
sity in favour of the Mo model, it makes more sense to ex- 
press the field strength estimation in terms of the upper 
limit of the credible regions. The null profiles are a good 
representation of what our data would look like in the ab- 
sence of a stellar magnetic field. For the whole N dataset, 
the upper limit of the 95.4% credible region is 144 G, well 
below the inferred dipole strength [r^ 1 kG) from the Stokes 
V observations. For the 99.0% credible region, we can say 
that the probability that an undetected field would have a 
dipole strength of more than 300 G is only 1%. The noise 
scaling does not significantly change the credible regions. 



Interestingly, the dipole strength inferred from the low s/n 
Stokes V observations is around 600 G, and the odds ratio 
favours the magnetic model. From the low s/n null N profile 
observation, we can see that a field with a dipole strength of 
the order of the 95.4% credible region upper limit (502 G) 
can indeed be detected. 



6 CONCLUSION 

In this paper, we have described a method based on Bayesian 
statistics to infer the magnetic properties of stars observed 
spectropolarimetrically in the context of large surveys like 
the Magnetism in Massive Stars project. This approach is 
well-suited for stars for which the stellar rotation period, 
and therefore the rotational phases of the small number of 
observations, are not known. 

The model used to predict the expected Stokes V 
profiles is that of an oblique dipolar magnetic field, 
parametrised by the field strength at the pole, the incli- 
nation of the rotational axis, the obliquity of the magnetic 
axis with respect to the rotational axis and the rotational 
phase (which is allowed to take any value for each individual 
observation). In the present case, the calculations are per- 
formed under the weak- field approximation, although any 
polarised spectral synthesis code can in principle be used 
with the Bayesian algorithm. 

The result of the analysis is a multidimensional poste- 
rior probability density that describes the relative likelihood 
of models spanning the parameter space of the dipolar field 
model. We have used synthetic observations to explore the 
behaviour of the Bayesian algorithm under ideal and realistic 
conditions. In the case of an ideal non-detection, the poste- 
rior probability density for the field strength has the form 
of a decreasing exponential, the extended tail of the field 
strength distribution being due to a specific family of dipole 
orientations. This tail becomes less extended with the addi- 
tion of multiple observations, as some of these specific dipole 
orientations only present low-amplitude or null Stokes V 
profiles over a restricted range of phases. However, the pos- 
sibility that the dipole is aligned with the rotational axis and 
seen equator-on {i = 90° and /3 = 0°) always remains as this 
configuration never produces any circular polarisation at any 
phase. When a detectable field is present, the probability 
distribution shows a sharp lower limit on the dipolar field 
strength, and a slow decrease towards higher field strengths, 
producing an asymmetrical distribution. With only one ob- 



Stellar magnetic parameters from Bayesian analysis 19 



servation, not much can be inferred about the dipole orien- 
tation. The longitudinal field indicates which pole is located 
on the visible hemisphere and the shape of the Stokes V 
profile indicates when some obliquity is necessary because 
the pole is located at a non-null rotational radial velocity. 

A particularly useful quantity that can be computed 
from the posterior probabilities is the so-called "odds ra- 
tio" , which compares the relative compatibility between the 
observations and the magnetic dipole hypothesis versus the 
non-magnetic hypothesis. By adding a magnetic signal cor- 
responding to the upper limit of the credible regions found 
for the ideal non-detection, we have explored the detection 
capability of the odds ratios. We find that fields correspond- 
ing to the 68.3% and 95.4% region are below detection, 
whereas those corresponding to the 99.0% and 99.7% re- 
gions are well detected with the odds ratios. In contrast, 
traditional diagnostics (detection probability and global lon- 
gitudinal field significance) only detect fields corresponding 
to the 99.7% region. 

By using a set of five realistic simulated observations, 
we have also shown that in the case of noise emulating the 
shape of a weak magnetic signal, it is possible to use the 
odds ratios to distinguish between noise and real signal by 
obtaining a small number of additional observations. Com- 
bining all the observations together, it is therefore possible 
to detect a weak magnetic field (with a strength correspond- 
ing to roughly the upper limit of the 95.4% credible region 
of the noise-only case) under the detection capabilities of 
the traditional diagnostics. We have therefore shown that 
the odds ratio is an powerful quantitative indicator of which 
undetected stars in a survey should be re-observed in prior- 
ity. 

In most applications where a field is indeed detected, the 
resulting probability density will generally be marginalised 
for the dipole strength, as only limited information can be 
obtained for the rotational inclination and the magnetic 
obliquity from a few observations (the magnetic geometry 
is recovered by the most probable inclination and obliquity, 
but the credible regions are quite extended) . We have shown 
that the dipole strength probability distribution provides a 
reasonable estimate of the field strength. 

We have applied our method to real spectropolarimet- 
ric observations of the magnetic B-type star LP Ori. The 
dataset consists of 3 ESPaDOnS and 4 Narval observations, 
of various signal-to-noise ratios. A magnetic signal is indeed 
detected by the odds ratios, even when only considering the 
low s/n observations where the traditional diagnostics do not 
detect the magnetic field. Using all the available spectra, we 
used the median of the marginalised posterior probability 
density, as well as the 68.3% credible region, to infer a dipo- 
lar field strength of 911 1244 G- Although the probability 
density for the obliquity and inclination do not provide any 
tight constraints, geometries for which the negative mag- 
netic pole spends little or no time on the visible hemisphere 
are preferred, since the dataset consists of positive or nearly 
null longitudinal field measurements. We also performed our 
analysis on the diagnostic null spectra, which resulted in a 
non-detection. The null profiles provide an useful verifica- 
tion of spurious signals and also provide an estimate of the 
field strengths that would be detectable with data of such 
quality. For example, when considering only the low s/n ob- 
servations, the Stokes V profile analysis yielded the detec- 



tion of the 600 G dipolar field and the diagnostic null 
profile analysis yielded an upper limit of ~ 500 G for the 
95.4% credible region. 
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Table 11. Line mask used for the LSD profile construction of LP Ori. The given values are the central wavelength, the ion, the unbroaden 
line depth with respect to the continuum and the effective Lande factor. 
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